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Abstract 

Symmetry breaking across phase transitions often causes changes in selection rules and emergence 
of optical modes which can be detected via spectroscopic techniques or generated coherently in 
pump-probe experiments. In second-order or weakly first-order transitions, fluctuations of the 
order parameter are present above the ordering temperature, giving rise to intriguing precursor 
phenomena, such as critical opalescence. Here, we demonstrate that in magnetite (FeaCH) light 
excitation couples to the critical fluctuations of the charge order and coherently generates structural 
modes of the ordered phase above the critical temperature of the Verwey transition. Our findings 
are obtained by detecting coherent oscillations of the optical constants through ultrafast broadband 
spectroscopy and analyzing their dependence on temperature. To unveil the coupling between the 
structural modes and the electronic excitations, at the origin of the Verwey transition, we combine 
our results from pump-probe experiments with spontaneous Raman scattering data and theoretical 
calculations of both the phonon dispersion curves and the optical constants. Our methodology 
represents an effective tool to study the real-time dynamics of critical fluctuations across phase 
transitions. 

* Corresponding author: fabrizio.carbone@epfl.ch. 
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INTRODUCTION 


Magnetite (FesC^) is the earliest discovered magnetic material. Nowadays, it is widely 
used in recording devices and as a catalyst [1, 2], and its fast reaction to light excitation 
holds promise for future applications in oxide electronics [3, 4], At room temperature, the 
ground state of magnetite is half-metallic ferrimagnetic and the crystalline structure is an 
inverse cubic spinel, in which tetrahedrally coordinated A-sites are occupied only by Fe 3+ 
ions, while octahedrally coordinated B-sites are occupied by an equal number of nominal 
valence Fe 3+ and Fe 2+ ions. At the Verwey temperature, Ty (= 116 K in our sample), 
important discontinuous structural and electronic modifications take place. In particular, 
the crystal symmetry lowers from cubic Fd3m to monoclinic Cc and the dc conductivity 
decreases by two orders of magnitude. The insulating state is the result of the long-range 
charge and orbital order (CO-OO) that emerges among the B-sites in the low-temperature 
(LT) phase. 

Originally, Verwey and Anderson postulated a purely electronic mechanism for the trans¬ 
formation [5, 6], which implies an alternate occupation of the a/4-spaced (001) c lattice planes 
(where a is the lattice parameter) by Fe 2+ and Fe 3+ ions below Ty, and a random distribu¬ 
tion of the latter cations above Ty, allowing for electronic transport. Although the above 
model proves correct to a first approximation, recently a far more tangled charge and orbital 
arrangement was refined below Ty, comprised of excess electrons delocalized over chains of 
three B-type Fe ions, termed trimerons [7-9]. The trimerons are examples of molecular 
polarons in which the electronic distribution is interrelated with the shortening of the Fe- 
Fe distances. Interestingly, extensive research on precursors of the phase transition (PT) 
substantiates the persistence of local dynamical CO-OO in the form of fluctuating trimeron 
complexes from Ty to room temperature [10-17]. However, the nature of the short-range 
order in the pre-transition region and the relation between the critical correlations which 
appear gradually above Ty and the discontinuous modifications at Ty remain contentious. 
Sample dependent effects on the critical phenomena prevent more general considerations 
[18, 19]. 

Optical probes, such as reflectivity [20], ellipsometry [4], infrared spectroscopy and spon¬ 
taneous Raman scattering [21, 22] have disclosed important information on the electronic 
and structural properties of magnetite relevant to the mechanism of the Verwey transition. 
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Transfers of spectral weight among different energy regions of the optical spectrum are sig¬ 
natures of the strongly correlated nature of magnetite [23], which manifest themselves both 
in the discontinuous change of the electronic structure at Ty and in the precursor variations 
which progressively preempt the Verwey transition. The absorption structures below 2.5 eV 
were generally attributed to interband transitions between Fe 3 d states, although alternative 
interpretations were proposed, at variance on whether only B- or also A-sites are involved 
[4, 20, 24, 25]. The appearance of a rich spectrum of additional infrared- and Raman-active 
modes below Ty reflects the symmetry lowering of the crystal structure which characterizes 
the Verwey transition [21, 22], 

The structural transformation of magnetite is a complex process, in which small atomic 
displacements from the prototype to the distorted structure are the result of the concomitant 
condensation of multiple phonons of comparable importance, rather than dominant contri¬ 
butions from few critical modes. The details of the distortions accompanying the Verwey 
transition have been elucidated by means of X-ray and neutron scattering [7, 9, 26-28]. In 
particular, below Ty, subsets of modes at the T-, A-, X- and W-point are frozen into the 
static deformations of the low-symmetry structure. Verwey and Anderson’s scenario, and 
subsequent Cullen and Callen’s theory [29] only account for electron-electron interactions 
and thus cannot explain the structural transformation. The involvement of the structural 
degree of freedom in the transition process was first put forward by Yamada, who pro¬ 
posed that the condensation of fluctuations of the charge density coupled to a A 5 symmetry 
phonon is the main driving force for the metal-insulator transition (MIT) [30]. More recent 
first principle computations suggest that both the A 5 and X 3 modes are the primary order 
parameters of the structural transformation, and the strong coupling of the X 3 mode with 
the conduction band electrons primarily causes the insulating state below T v [31, 32], 

In Yamada’s model, the modes of the charge density are treated as classical Ising vari¬ 
ables and the transformation as an order-disorder process. Therefore, the spectrum of the 
electronic excitations of the order parameter is essentially a relaxation response centered at 
zero frequency, whose width is determined by the inverse relaxation time. A radically differ¬ 
ent model considers a band of quasi-one-dimensional fermionic quasiparticles, generated by 
the correlations in the system, and interprets the PT in terms of a Peierls instability [33]. 
In this case, one expects a richer spectrum of electronic excitations, namely, quasiparticle 
transitions across the Peierls gap and a collective mode arising from amplitude fluctuations 
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of the order parameter at finite frequency [34], 

To map the dynamical evolution of the structural modes across the Verwey transition 
and their interplay with the electronic excitations, we performed out-of-equilibrium optical 
spectroscopy experiments. We found that the critical fluctuations of the charge-density- 
wave (CDW) couple to specific ionic motions, giving rise to precursor effects above the 
critical temperature. In particular, in this paper we show that an optical pump pulse can 
transiently promote time-dependent electron-phonon interactions in a sample of natural 
magnetite and generate coherent phonons which, by symmetry, should be forbidden above 
the charge ordering temperature. The launch of these forbidden phonons needs the assistance 
of fluctuations of the ordered phase. This is different from previously proposed mechanisms 
of coherent phonon generation [35, 36] in which a time-dependent force acts on the ions 
without the assistance of other low-lying modes. 

Ab initio calculations of the phonon dispersion curves in the monoclinic structure allow 
us to identify the collective modes that are impulsively photoexcited with the X 3 (13.9- 
15.6 meV), A 2 (18.5-19.6 meV) and T 2 g mode (25.9 meV). Our assignment of the coherent 
phonons is further complemented by theoretical predictions of the Raman matrix elements 
(RME). The correspondence between data and computations serves as an anchor point to 
rediscuss the attribution of the optical absorption peaks at low energy. As a by-product, 
we clarify the origin of the transitions at optical frequencies debated in literature. Further¬ 
more, our results provide a measure of the amplitude and persistence time of the critical 
fluctuations above Ty. 

EXPERIMENTS 

A natural single crystal of magnetite was cut and polished along the (110) c plane. In the 
Supplemental Material (SM), we report a thorough characterization by means of resistivity 
and ac magnetic susceptibility experiments, revealing the typical fingerprints of high-quality 
single crystals. 

Figure 1 displays the steady-state optical conductivity of magnetite in the low energy 
region, where temperature dependent changes of the spectral weight take place. As Ty 
is crossed, an optical gap opens below 0.2 eV, signaling the onset of the insulating phase. 
Three main features are respectively centered at 0.6, 2 and 5 eV. Based on density functional 
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theory (DFT) calculations [24], Park and co-workers first attributed such features to intersite 
transitions between B-type Fe t 2 g states, between A- and B-sites and the charge-transfer 
(CT) between O 2 p and Fe 3 d states, respectively [20]. Advanced DFT+U calculations 
revised the latter assignment, indicating that the absorption feature at 2 eV does not involve 
A-type Fe sites, but only B-type Fe sites with different valence states [25]. Nevertheless, no 
meaningful intensity was predicted for this peak, which makes the assignment put forth in 
Ref. 25 open for discussion. 

In our pump-probe experiments, the photoexcitation energy was tuned to the CT region 
(3.1 eV) and the probe band spanned the range of the d-d transitions (1.7-2.9 eV). The 
sample was fixed to the copper cold-finger of a closed-cycle He cryostat, with a base pressure 
of 10 -8 mbar. The output of an amplified TkSapphire laser delivering 50 fs pulses centered 
at 1.55 eV at a repetition rate of 6 kHz was split into a pump and a probe beam. The pump 
beam was used to photoexcite the material, after doubling the fundamental frequency in 
a beta-barium borate crystal, whereas the probe beam was focused on a CaF 2 crystal to 
generate broadband pulses from 1.7 to 2.9 eV. The reflectivity spectra were recorded by a 
fiber-coupled spectrometer. A chopper modulated the train of pump pulses at 3 kHz, so that 
at every delay time the difference between the pumped and unpumped reflectivity spectrum 
was measured, yielding a sensitivity down to 10~ 4 . Further details of the set-up are described 
in Refs. 37 and 38. The delay time between the two pulses was varied to obtain a sequence 
of time points representative of the photoinduced dynamics of the material. Spontaneous 
Raman scattering at 2.4 and 3.1 eV excitation energies was also carried out to discriminate 
between equilibrium and non-equilibrium effects (see SM for further details). 


RESULTS 

Ultrafast broadband spectroscopy 

The spectrum of the transient reflectivity as a function of pump-probe delay time, 
AR/R (E, t), at 10 K and ~1 mJ/cm 2 fluence is shown in Fig. 2a. The AR/R tempo¬ 
ral traces averaged over 0.2 eV wide energy intervals are displayed in Fig. 2b. Damped 
coherent oscillations lasting over 1.5 ps superimposed to a multi-exponential decay clearly 
emerge in all the time traces of Fig. 2b and in the colormap of Fig. 2a. To analyze these 
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oscillations, the temporal trace averaged between 2.0 and 2.3 eV probe photon energy is 
fitted with the sum of multi-exponential decays and two distinct damped coherent oscilla¬ 
tions, convolved with a Gaussian response function. The data and the fitting functions are 
displayed in Fig. 2c. The beating between the two oscillations shown in Fig. 2d is thereby 
isolated. 

The same fit analysis is repeated at different temperatures and the same fluence. The 
temperature dependence of all oscillations is displayed in Fig. 3a-d. Figure 3e-h shows the 
evolution of the oscillation amplitude, as obtained by Fourier transforming the time traces. 
According to the fit results, the energies of the two main modes are respectively in the 
13.9-15.6 meV and 18.5-19.6 meV ranges. 

As illustrated in Fig. 3b and 3d, in the cubic phase, AR/R also shows a weak coherent 
oscillation corresponding to an energy of 25.9 meV, visible in a different spectral range 
(< 1.9 eV). Remarkably, although the intensity of the 13.9-15.6 meV and 18.5-19.6 meV 
oscillations decreases abruptly above Ty, it does not vanish up to 140 K (see Fig. 3e and 
3g). Instead, the 25.9 meV oscillation is only visible in the high-temperature (HT) phase 
and progressively weakens with increasing temperature (see Fig. 3f and 3h). 

Next, in Fig. 4a-d, we show the dependence of AR/R on the fluence at 120 K, right 
above the transition temperature. As evidenced in Fig. 4d, the 13.9-15.6 meV and 18.5-19.6 
meV oscillations are completely quenched with increasing fluence above ~1 mJ/cm 2 . In 
contrast, the intensity of the 25.9 meV oscillation increases with fluence and saturates above 
~4 mJ/cm 2 . 


Mode assignment and phonon calculations in P2/c symmetry 

To assign these modes, a direct comparison of the FT of the coherent oscillations and 
the spontaneous Raman spectra with the theoretical predictions is carried out in Fig. 5a- 
d. The phonon dispersion curves obtained from ab initio calculations in the cubic phase 
below 27 meV are represented by red lines in Fig. 5c. Figure 5d shows the spontaneous 
Raman spectra at room temperature for 2.4 and 3.1 eV excitation energies, and the FT of 
the coherent oscillations detected in our pump-probe experiments. The only Raman-active 
phonon at the zone center in the relevant energy region is the T 2 g mode. 

In the monoclinic phase, the structural PT causes a quadrupling of the unit cell, resulting 


7 



in new phonon branches, as shown by the blue lines in Fig. 5b. As a consequence, various 
modes are folded to the zone center and become potentially observable by visible-light Raman 
scattering, as illustrated by the horizontal lines in Fig. 5a. Figure 5a also displays the 
spontaneous Raman spectra of our single crystal at 5 K for 2.4 and 3.1 eV excitation energies, 
and their time-resolved counterpart at 10 K. The energy ranges of the coherent oscillations 
as inferred from the fit analysis are highlighted by shaded areas. 

The phonon dispersion curves were calculated in the monoclinic structure (with the ap¬ 
proximate symmetry P2/c) using the ab initio direct method [39, 40]. The electronic struc¬ 
ture and atomic positions were optimized using the projector augmented-wave [41] and gen¬ 
eralized gradient approximation (GGA) [42] implemented in the VASP program [43]. The 
strong electron interactions in the Fe(3d) states were included within the GGA+U method 
[44], The Hellmann-Feynman forces were calculated by displacing all non-equivalent atoms 
from their equilibrium positions and the force-constant matrix elements were obtained in the 
112-atom supercell. The phonon dispersions along the high-symmetry directions in the re¬ 
ciprocal space were calculated by the diagonalization of the dynamical matrix. The phonon 
dispersions in the HT cubic phase were studied using the same approach and discussed in 
Refs. 31 and 32. 

We observed that for crossed pump and probe polarizations the signs of the oscillations 
visible in the monoclinic phase do not change, in contrast to the behavior typical of B g 
symmetry phonons. Therefore, we rule out the LT B g modes from the potential candidates 
for the assignment of these oscillations. On the basis of the comparison presented in Fig. 5, 
the 13.9-15.6 meV and 18.5-19.6 meV oscillations are ascribed to the LT A g modes that 
are closest in energy, which originate from respectively the lowest-energy X 3 and A 2 modes 
in the HT phase. Instead, the 25.9 meV oscillation observed above Ty is unambiguously 
identified by its energy with the T 2g mode of the cubic structure. 

Let us note that the peaks in the FT of the coherent oscillations extend in a broad 
energy range from 11.5 to 21.0 meV. The occurrence of oscillation dephasing and energy 
renomalization from anharmonic effects may explain the large width of such features and 
the lower energy of the oscillations compared to their spontaneous Raman counterparts. 
Our calculations in the approximate symmetry P2/c do not fully account for the effects of 
C0-00 on the phonon energies, which are thus underestimated, incidentally resulting in a 
better agreement with the time-resolved data than with the spontaneous Raman spectra. 



The cubic A 5 mode gives rise to a monoclinic A g mode in the energy range intermediate 
between the 13.9-15.6 meV and 18.5-19.6 meV oscillations, which may cast donbts on the 
above assignment based on energy considerations. Such doubts will be dispelled in the 
following. 


RME analysis and optical constant calculations in P2/c symmetry 

To corroborate our assignments of the oscillations, we benefit from the broadband nature 
of our optical probe and compare the energy dependence of the experimental RMEs of the 
oscillations to those obtained from our ab initio calculations. Singular value decomposition 
(SVD) of the AR/R matrix was performed to separate the energy dependence of the incoher¬ 
ent relaxation of the excitations from that of the collective coherent modes corresponding to 
the oscillations. Qualitative information on the coupling between structural distortions and 
electronic excitations was thereby gained. Our algorithm decomposed the AR/R matrix into 
the sum of outer products between orthonormal spectral and temporal vectors, also termed 
canonical traces, weighted by singular values. The sum was truncated to account only for 
the physically relevant terms, thus reducing the noise level. The canonical traces were fitted 
with model functions (see Supplementary Fig. 4). By that, as shown in Fig. 6 a and 6 b, the 
incoherent and the coherent contribution to AR/R were reconstructed. The same analysis 
combined with a fit using a Lorentz model delivered the energy profiles of the differential 
permittivity associated with the two oscillations at the maximum displacement from the 
equilibrium positions, depicted in Fig. 6 c-e and in Fig. 7b, in a wider energy range. A 
thorough description and application of our SVD-based method is available in Ref. 37. 

To obtain the theoretical RMEs, we followed the procedure described in Ref. 37. The 
frequency-dependent dielectric function was calculated for a series of monoclinic unit cells 
(with approximate symmetry P2/c) whose atoms were displaced according to the eigen¬ 
vectors of the phonon modes in the same 11.5-21.0 meV energy range as the oscillations 
observed in our experiments. These are the A g modes of the monoclinic structure that orig¬ 
inate from the cubic X 3 , A 5 and A 2 counterparts. We computed the RMEs as the finite 
differences of the dielectric function associated with the structural distortions in the linear 
displacement regime. The self-consistent charge density of the electronic ground state was 
obtained from DFT calculations in the GGA+U approximation implemented in the QUAN- 
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tum ESPRESSO package [44, 45]. Subsequently, the frequency-dependent dielectric function 
was calculated with the linear response method within the random phase approximation 
using the yambo package [46] (see the SM for computational details). For values of the 
on-site Coulomb repulsion U and the Hund’s exchange coupling J of respectively 4 eV and 1 
eV, there is qualitative similarity between the computed and experimental RMEs for the X 3 
and A 2 modes, whereas they are in complete disagreement for the A 5 mode (see Fig. 6 c-e). 
Notice that the computed RMEs depend on the energy locations of the transitions which, 
in general, are only approximately predicted by ab initio calculations. This may explain 
the relative energy shift of data and predictions for the X 3 mode (see Fig. 6 c), but not the 
inconsistency of the theoretical RME for the A 5 mode with both measured RMEs (see Fig. 
6 d). Thus, our RME analysis provides an anchor point to our assignment of the oscillations 
in the foregoing following from energy considerations alone. 

In addition, as illustrated in Fig. 7a, our calculations qualitatively reproduce the char¬ 
acteristic features of the optical conductivity at 15 K, except for the position of the lowest- 
energy peak. We note that the agreement with the experimental data did not improve for 
other possible choices of the U parameter (see Supplementary Fig. 5). Considering the 
reasonable consistency of our predictions to the experimental RMEs in the probed energy 
range, we re-examine the absorption structure of magnetite in the LT phase to better sub¬ 
stantiate the assignment of the interband transitions. As illustrated in Fig. 7a, the total 
spectrum was decomposed into the main contributions from different interband transitions 
(our complete analysis is presented in Supplementary Fig. 6 ). Individual electronic tran¬ 
sitions were discriminated based on the symmetry of the energy bands corresponding to 
different ions, summarized in Fig. 7c. Three separate components dominate the debated 
feature around 2 eV (represented with different colors in Fig. 7a), namely, the CT between 

(i) the minority-spin t 2 g and e g states of respectively nominal valence Fe 2+ and Fe 3+ B-sites, 

(ii) the minority-spin t 2 g states of the same ions and (iii) the top of the highest filled band 
and the bottom of the lowest empty band formed by respectively B- and A-type Fe ions in 
the majority-spin channel. Earlier calculations [25] attributed the peak centered around 2 
eV only to the first of the above three excitations purely based on energy considerations, 
albeit with negligible calculated intensity. In contrast, our assignment relies on the com¬ 
puted spectral weights, which reveal the additional involvement of the A-sites. Furthermore, 
we found that the excitation of the t 2g electrons across the bandgap extends in the energy 
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region relevant to the 2 eV feature in the form of a satellite peak. 


CDW-assisted excitation of the LT modes above Ty 

In the HT phase, the X 3 and A 2 modes are not at the zone center, thus they should not be 
accessible by low-momentum-transfer experiments such as Raman scattering (see Fig. 5c). 
Here, we show that the forbidden modes can become active with the assistance of CDW 
fluctuations. Such precursor effects can be expected in PTs, except for transitions which are 
mean-field-like or have a strong first-order character. In magnetite, such fluctuations have 
been demonstrated to persist above T v [10, 17] and provide the (0,0,1) and (0,0,|) wave- 
vector components (in units of 2tt/ a) for momentum conservation, in the Raman process 
pictorially represented in Fig. 8. 

We model the light pulses as a time-dependent electric held E. with an oscillation in the 
range of visible light frequencies, corresponding to the laser central frequency c Ol, and a 
slowly-varying modulation £, representing the pulse shape, E(f) = £(t)e lU}Lt . We introduce 
boson creation (a qA ) and destruction (a q A) operators for phonons of branch A and momentum 
q, and the corresponding canonical coordinate, Q^x = (h/2u cl \) 1 ^ 2 (a q a + at A ). The order 
parameter can be either a bosonic or a classical variable, depending on whether one considers 
the PT as a Peierls instability, with quantum fluctuation effects [47], or an order-disorder 
process, described by a 0 4 theory, with a central peak-like dynamical response [48]. 

The time-dependent perturbation introduced by the pump electric held can be expressed 
by the following two-mode Raman Hamiltonian, which couples CDW fluctuations and a 
phonon mode, thus acting as a time-dependent electron-phonon interaction controlled by 
the photoexcitation, 


Hr ^ ^ 9q\(t)pgQ—qX 

( 1 ) 

qA 



( 2 ) 


Here, p q is the canonical coordinate for the charge fluctuations with momentum q, V is the 
sample volume, and d 2 x/dp^dQ-^x is the Raman tensor, expressed as the second derivative 
of the dielectric susceptibility evaluated at the pump frequency, lol■ Such Raman tensor 
is a straightforward generalization of those in Refs. 35, 38, 49. For simplicity, we assume 
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that bare electron-phonon interactions are already incorporated in the definitions of the 
fluctuation operators, which are quasiparticle operators. For the charge fluctuations, the 
derivative can be evaluated with the use of a Lagrange multiplier, as in Ref. [38]. 

The response produced by the above interaction can be computed as follows, using the 
Kubo formula, in the response regime linear in fluence, 

(PqQ-qA>(t) = J df'([pq(f)Q_q A (£),<2qA(0P—q(O])0q(O 

= f dt'R(t — t')gq(t'). (3) 

J — oo 

where (...) indicate thermal and quantum averages. If charge fluctuations behave as a 
classical order parameter, p q can be taken out of the commutator, yielding, 


(PqQ-q)W — 

rt 


( 4 ) 


... / \ . A ,e ^ t '^ Tph sinlujait — t')] , A 

dt (p q (t)p_ q (f ))--0q(* ) 




where r p /, is a phenomenological phonon lifetime and o; q is the phonon frequency. For 
notational simplicity, we drop the phonon branch index A and consider a single phonon. 

An impulsive form of p q (f) produces sine-like oscillations of the charge-ion correlation, 
with an initial strength given by the FT of the squared amplitude of the order parameter 
fluctuations, (p q (0)p_ q (0)), and an envelope, (p q (f)p_ q (f')), which carries information on 
their time correlation. In the case of a step form of g q (f), the oscillation becomes cosine-like. 
Our description represents a generalization of the mechanism proposed in Ref. 36 which can 
be referred to as fluctuation-assisted stimulated Raman scattering (FASRS). 

The launch of the excitation modulates the optical properties through their dependence 
on the dielectric susceptibility, 


$x(u,t) = Y QpQQ ^)M-q)W- ( 5 ) 

To apply the FASRS mechanism to the present case, we assume a CDW instability 
occurring at the X-point of the Brillouin zone (BZ), corresponding to a charge stacking along 
the z direction, described by cos(27rz/a), which yields a charge occupancy (with respect to 
the average) —5, 0, S, 0 for the B-type Fe planes separated by a/4. One can additionally 
consider CDW fluctuations with twice the periodicity, which can assist the launch of modes 
at the A-point of the BZ, and can be described with the same method. 
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The dominant CO in real space, with Q=(0,0,27r/a) wave-vector, corresponding to the 
X-point of the BZ, can be described as, 

Sp(r) = </>i(r) cos(Q.r) + 0 2 (r) sin(Q.r). 


where the trigonometric factors account for the rapidly varying part of the CO, whereas 
(0i, 02) is a smoothly varying real order parameter. Expanding the energy in a power series 
of the order parameter and of its gradients yields a coarse-grained Landau energy functional, 
describing the PT [48]. 

Close to the ordering temperature, critical fluctuations of the density appear strongly 
peaked at the critical wave-vector and we can approximate the smooth momentum dependent 
quantities in Eq. (5) to their value at the ordering wave-vector. The response reads, 

d 2 x , > N 2 


&x(w,t) = 


dp Q 8Q. Q M 2 X 


( 6 ) 


, / / , / \ i , A . e tn >/ T v h sinLWf — t')] . 

dt V (0aO (t) 0aO (t ) ) - ~9Q (t ) 

a=l,2 U <* 


with N the number of sites and (0 a o(^)0oo(0)) the order parameter autocorrelation at r = 0. 
Neglecting the fluctuations in Landau theory, the autocorrelation is proportional to the 
square of the order parameter, which increases linearly in Ty-T with decreasing tempera¬ 
ture, as schematically illustrated in Fig. 9 (blue line). This corresponds to the notion that 
the modes become active in the broken symmetry phase upon folding from the boundary 
to the center of the BZ. If the fluctuations are taken into account, the results are differ¬ 
ent. The equal-time correlation determines the initial amplitude of the oscillations. If one 
computes the correlation at the Gaussian level, one obtains the red line in Fig. 9, showing 
that the fluctuations make the modes active even above Ty. Furthermore, if the oscillations 
are shorter than the phonon lifetime, it means that the decay is dominated by the factor 
(0ao (t)0 ao (^)) Eq. (6) and one obtains information on the dynamics of the order param¬ 
eter. For a weakly first-order transition, as in the present case, we expect a small jump of 
the intensity at the critical point, but a qualitatively similar behavior. 


DISCUSSION 

Over the last decades, revived attention has been focused on the precursors of the PT 
which appear far above the critical temperature. For instance, disk-like diffuse scattering 
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anticipates the Verwey transition starting from room temperature, as a result of dynamical 
correlations among trimeron complexes, sustained in the cubic phase [11, 12, 17]. Recently, a 
lattice dynamical study found anomalies of transverse-acoustic phonons at incommensurate 
wave-vectors, which originate from their strong coupling to the same fluctuations of the 
charge density [16]. In general, the experimental results provide extensive evidence to the 
thermal population of degenerate excited state of short-range order prior to the Verwey 
transition. Degeneracy originates from the equivalent orientations of the order pattern 
in cubic symmetry, similar to crystal microtwinning in the monoclinic phase, but on a 
nanometer scale and in dynamical conditions. 

Within the same framework, the excitation and detection of two modes of the monoclinic 
structure above the critical temperature in our pump-probe experiments can be rationalized 
considering the effects of critical fluctuations and provide information on their dynamics. 
Eqs. (4) and (6) express the fact that, in the presence of a long-lived CDW fluctuation, a 
lattice vibration with finite momentum can be excited by light, with a coherence time limited 
by the persistence time of the CDW fluctuation. Practically, the frequency of the modes 
does not change across the critical region. This means that the modes are rendered active by 
classical fluctuations of the order parameter, as expected in Yamada’s model, consistent with 
the central peak observed in inelastic neutron scattering [10-12, 17]. Instead, if a Peierls 
mechanism was playing a significant role, one would expect that the electronic spectrum of 
the modes affects the response. For example, if one approximates the electronic spectrum 
to a single pole at cjcdw, the oscillations should appear at cuq ±cjcdw, with a strong energy 
dependence on temperature, which we do not observe. 

By means of our SVD algorithm, we were able to determine the experimental profiles of 
the R.MEs as a function of probe photon energy. A comparison with ab initio calculations 
allowed us to unambiguously identify the modes. To the best of our knowledge, this assign¬ 
ment methodology is introduced here for the first time and has the advantage to rely on the 
phonon eigenvector instead of the phonon energy, thus providing more precise and specific 
information than other optical approaches. 

Our time-resolved experiments directly access the decay time of the coherent oscillations. 
According to Eq. 6, this is limited by the phonon lifetime r p h or the persistence time of the 
CDW fluctuations, whatever is shorter. If the fluctuations of the order parameter are the 
limiting factor, we expect that the lineshapes of the modes in the energy domain become 
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broader above the critical temperature, as it is indeed the case (see Fig. 3e). Prelimi¬ 
nary results from our inelastic neutron scattering [50] and spontaneous Raman scattering 
experiments (see Supplementary Fig. 3) indicate that the phonon lifetimes (> 600 fs) are 
significantly larger than the time constants of the coherent oscillations in accord with this 
result. Therefore, our measurements provide an estimate of the relaxation time of the fluc¬ 
tuations of the order parameter in the 130-140 K pre-transition region, tcdw ~ 240 fs. 

Gaussian fluctuations result in a cusp anomaly of the oscillation amplitude at the critical 
temperature (see Fig. 9). Unfortunately, our temperature and energy resolution are not 
high enough to test this hypothesis. In any case, the persistence of fluctuations above the 
ordering temperature (see Fig. 3g) is in qualitative agreement with the prediction of the 
fluctuation model. The intensity drop with increasing fluence for the A 2 mode (see Fig. 4d) 
is probably due to temperature detuning from the critical region, and represents a hint of 
the fast decrease in intensity above the ordering temperature expected from the fluctuation 
model. 

In contrast to our non-equilibrium experiments, in our spontaneous Raman measurements 
at 2.4 and 3.1 eV excitation energies, no forbidden mode was detected below 27 meV in the 
disordered phase (see Supplementary Fig. 3). This is indicative of the different sensitivity 
of the two techniques to low-energy modes, but could also suggest that non-linear effects 
enhance the coherent response in our pump-probe experiments. One can easily check that the 
time-dependent electron-phonon interaction promotes a transient stabilization of the ordered 
phase, which is equivalent to an increase of the effective ordering temperature in steady- 
state conditions. Possibly, the above effect amplifies the coherent response, by inducing a 
breaking of symmetry limited in time and confined in space. Experimentally, to validate 
this hypothesis, one would need to carry out a thorough study of the fluence dependence 
of the coherent response in the proximity of the critical temperature. Such an experiment 
requires accurate temperature control across the specific heat anomaly in the critical region 
and high enough signal-to-noise ratio to detect small oscillations at low fluences. 


CONCLUSIONS 

Typically, ultrashort laser pulses are used to melt the order parameter of the low- 
symmetry phase and thus establish the electronic and structural properties of the disordered 
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phase [3, 4]. Here, we discussed a more subtle effect, with implications on the coherent con¬ 
trol of transition-metal oxides. In particular, we demonstrated that impulsive photoexcita¬ 
tion at 3.1 eV, above the energy onset for O 2p - Fe 3 d CT transitions, transiently promotes 
the coupling between lattice vibrations at finite momentum and fluctuations of the electronic 
order parameter above the ordering temperature of the Verwey transition in magnetite. In 
particular circumstances, the optical control of the electron-phonon interactions may even 
induce order. 

Fluctuations of the order parameter are not easy to measure and often require neutron 
scattering techniques, generally limited in energy resolution. We developed a method able 
to probe the dependence of the amplitude and persistence time of critical fluctuations on 
temperature and light excitation. Our method demonstrated the capability to gain insight 
into the critical softening of fluctuations in correspondence of PTs. Moreover, our approach 
enabled us to describe the atomic motions associated with the critical fluctuations and their 
interplay with the electronic degree of freedom. This further elucidated the microscopic 
mechanism of the MIT and the interpretation of the optical response of magnetite. 
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Figure 1 | Steady-state optical properties. Real part of the optical conductivity at 15 
and 300 K from the ellipsometry measurements of Ref. 4. Lines are dashed where curves are 
extrapolations based on fits to the data. The blue arrow and the shaded area respectively 
correspond to the energy of the impulsive photoexcitation and the probed range in our 
time-resolved experiments. 
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Figure 2 | Response to the impulsive photoexcitation in the monoclinic phase. 

(a) Color-coded map of differential reflectivity as a function of pump-probe delay time and 
probe photon energy at 10 K and ~1 mJ/cm 2 fluence. (b) AR/R time traces averaged over 
0.2 eV wide windows centered on the horizontal lines of panel a. (c) Fit analysis of AR/R 
integrated between 2.0 and 2.3 eV. Fitting functions are displayed in different colors, (d) 
Coherent response isolated through fit analysis of AR/R integrated between 2.0 and 2.3 eV. 
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Figure 3 | Temperature dependence of the coherent effects. (a,b) Differential 
reflectivity as a function of pump-probe delay time averaged over 2.0-2.3 eV and 1.72- 
1.92 eV energy ranges at different temperatures. (c,d) Oscillatory component singled out 
from the time traces of panel a and b by subtracting the non-oscillatory transient. (e,f) 
Power spectrum of the coherent response from a Fourier transform (FT) of the oscillatory 
component. Blue (Red) vertical lines indicate the energy of the A g (T 2 g ) phonon modes at 
the T-point in P2/c ( Fd3m ) symmetry from ab initio calculations. Shaded areas extend 
over the energy ranges of the beating modes in the LT phase as estimated based on the fit 
analysis. (g,h) Temperature dependence of the peak |FFT | 2 corresponding to the oscillations 
visible by eye. A vertical dashed line denotes Ty. 
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Figure 4 | Fluence dependence of the coherent effects at 120 K. (a) Differential 
reflectivity as a function of pump-probe delay time averaged in between 2.0 and 2.3 eV at 
different fluences. (b) Oscillatory component singled out from the time traces of panel a 
by subtracting the non-oscillatory transient, (c) Power spectrum of the coherent response 
from the FT of the oscillatory component. A red vertical line indicates the energy of the 
T 2g phonon mode at the T-point in Fd3m symmetry from ab initio calculations, (d) Tem¬ 
perature dependence of the peak |FFT| - corresponding to the oscillations visible by eye. 
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Figure 5 | Assignment of the phonon modes in pump-probe and spontaneous 
Raman experiments. (a,d) Comparison between our spontaneous Raman spectra at 
different excitation energies (colored lines) and |FFT| J of the coherent effects in our time- 
resolved experiments (black line), respectively in the monoclinic and cubic phase. Blue (Red) 
horizontal lines indicate the energy of the A g (T 2 g ) phonon modes at the T-point in P2/c 
(Fd3m) symmetry from ab initio calculations. Shaded areas extend over the energy ranges 
of the beating modes in the LT phase as estimated based on the fit analysis. (b,c) Ab initio 
calculations of the phonon dispersion curves, respectively in P2/c and Fd3m symmetry. 
Labels indicate the symmetry of the modes at the T-, A- and X-point in Fd3m symmetry, 
and the cubic counterparts of the Raman-active modes at the zone center in P2jc symmetry. 
Blue (Red) diamonds are referred to A g (B g ) modes. 
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Figure 6 | Time and energy dependence of the coherent response. Reconstruction 
of the incoherent (a) and the coherent contribution (b) to the AR/R matrix at 10 K and ~1 
nr J/cm 2 fluence through our SVD-based method, (c—e) Energy profiles of the contributions 
to the imaginary part of the differential permittivity, Ae 2 , from the oscillations at maximum 
displacement amplitude at 10 K and ~1 mJ/cm 2 (symbols), and calculated RMEs for the A g 
modes in the same energy region, in the monoclinic phase (lines). The cubic counterparts of 
the modes are specihed between parentheses. The computed RMEs correspond to maximum 
atomic displacements of about 5 x 10~ 2 A for the 15.3 meV and 18.8 meV modes and 1 x 10~ 2 
A for the 17.3 meV mode. 















Figure 7 | Ab initio calculations of the optical properties, (a) Comparison between 
the real part of the optical conductivity at 15 K from the cllipsometry measurements of Ref. 
4 and the theoretical predictions from our GGA+U calculations with U = 4 eV and J = 
1 eV. The main contributions from different interband transitions to the features around 
1.2 and 1.7-2.2 eV are plotted with colored lines and labeled according to the predominant 
character of the initial and final states, (b) Calculated RMEs for the modes corresponding 
to the oscillations in the low energy region below 3 eV. (c) Calculated partial density of 
states corresponding to inequivalent Fe and O ions in the P2/c subcell. The notation Bl-4 
refers to inequivalent positions in the P2/c subcell (see Ref. 7). The Fermi level (dashed 
line) is set to the top of the valence band. 
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Figure 8 | Two-mode photoexcitation mechanism. Schematic representation of the 
momentum conservation rule in the photoexcitation mechanism of the X 3 , A 2 and CDW 
fluctuation modes above Ty. The phonon dispersion curves in the cubic phase as measured 
by means of inelastic neutron scattering (symbols) [26] and predicted according to ab initio 
calculations (lines) are shown. The energy and momentum of the phonon and CDW fluctu¬ 
ations are represented by red and blue arrows, respectively. The combination of electronic 
and vibrational modes allows momentum conservation in the optical process. The shaded 
area extends over the energy range of the coherent oscillations in the cubic phase, in the 
vicinity of Ty, as estimated based on a FT analysis. 
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Figure 9 | Squared amplitude of the order parameter (fluctuations) as a function 
of temperature, below (above) the ordering temperature. The plotted curves are 
referred to a one-component </> 4 model at the Gaussian level. The results for a two-component 
model are qualitatively similar, but the details depend on the interactions among the modes. 
In these units, the intensity of the anomaly at Ty is regulated by the inverse Ginzburg length 

[48]. 
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51. SAMPLE CHARACTERIZATION 

The resistivity and ac magnetic susceptibility of our sample were measured as a function 
of temperature, respectively in the ranges 35-300 K and 5-150 K, to determine the Verwey 
temperature (Ty) and the transition order. The ac susceptibility measurements were per¬ 
formed using a commercial CryoBIND system. Resistivity and ac magnetic susceptibility 
data are displayed in Fig. SI. The resistivity data show that the transition is first order and 
Ty = 116 K, below typical values for highly stoichiometric samples, owing to the impurity 
content inherent to natural rocks of magnetite. As illustrated in Fig. Sib, the ac magnetic 
susceptibility of our sample exhibits features analogous to Refs. 51 and 52, referred to syn¬ 
thetic crystals. In particular, in accordance to the above studies, the magnetic response is 
independent of frequency in the critical region, whereas a pronounced frequency dependence 
characterizes the susceptibility anomaly at low temperature. The magnetization dynamics 
is ascribed to domain wall motion (DWM) according to Ref. 52. The drop in the real part of 
the ac magnetic susceptibility X’ originates from crystal microtwinning across the structural 
transition. In the monoclinic phase, the magnetic domains are coincident with the ferroclas- 
tic domains formed below Ty, which require higher fields to propagate via DWM, compared 
to their cubic counterparts. Therefore, the temperature for the discontinuity in X’ is under¬ 
stood as the critical temperature of the structural transition, and found to correspond to 
that of the metal-insulator transition. 

52. SPONTANEOUS RAMAN SCATTERING 

Figure S3 shows the spontaneous Raman spectra of magnetite as a function of temperature 
at 2.4 and 3.1 eV excitation energies. In the high temperature phase, the T 2 g mode is visible 
only for 2.4 eV excitation energy, which further substantiates the dependence on energy of 
its Raman matrix element (RME), as deduced from our time-resolved data. The energy 
of the T 2g mode is in agreement with Ref 21. The symmetry breaking across the Verwey 
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transition gives rise to a fine structure of Raman-active modes at the Brillouin zone center. 
The red arrows in Fig. S3 mark two peaks which appear below Ty and are not related 
to any peak in the high-temperature phase. On the basis of the analysis presented in the 
main text, the peak at 20.2 meV at 5 K is assigned to the A g mode that originates from 
the A 2 mode in the cubic phase. The dependence on temperature of its energy and width 
is depicted in Fig. S3c. A lower energy mode is barely visible at 15.6 meV at 5 K, only for 
2.4 eV excitation energy, owing to experimental limitations to the lower cut-off energy when 
using 3.1 eV light. Based on the same analysis, this mode is attributed to the lowest-energy 
X 3 mode in the cubic phase. 


S3. AB INITIO CALCULATIONS OF OPTICAL CONSTANTS AND RMES 

The yambo code for many-body calculations in solids [46] was used to compute the optical 
absorption spectrum within the linear response random phase approximation. The GGA+U 
ground state (GS) Kohn-Sham states calculated with the quantum espresso (QE) suite 
of programs [45] were used as input for the optical absorption calculations. The GS self- 
consistent field (SCF) charge density was obtained using norm-conserving pseudopotentials 
from the QE-pseudopotential library, a plane wave basis-set with 100 Ry energy cutoff and 
560 K-points sampling the first Brillouin zone. The experimentally observed charge ordering 
on the Fe B-sites was well reproduced with on-site Coulomb repulsion U = 4 eV and Hund’s 
exchange coupling J = 1 eV, yielding an insulating gap of 0.35 eV. As shown in Fig. S5, the 
optical constants obtained for U = 3.5 eV and U = 5 eV are in poorer agreement with the 
data compared to those calculated for U = 4 eV. 

The RMEs were computed for displacement amplitudes corresponding to once ( A = 1 ) 
and twice ( A = 2 ) the mode contributions to the structural distortions from Fd3rri to 
P2/c symmetry, starting from the equilibrium positions in the monoclinic structure. The 
maximum displacements of the atoms for A = 1 are about 8 x 10~ 3 A, 6 x 10~ 3 A and 
13 x 10 ” 3 A, respectively for the X 3 , A 5 and A 2 modes. The RME calculated for A = 2 and 
the X 3 (A 2 ) mode was multiplied by a factor 3 ( 2 ) for the comparison with the experimental 
RMEs (Fig. 6 c-e of the main text), which is equivalent to displacing the atoms by A = 6 
(4) if non-linear effects are negligible. Our results for A = 1 and A = 2 show non-negligible 
non-linear effects, which however remain small in the probed spectral range. 
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S4. SUPPLEMENTARY FIGURES 



FIG. SI. Resistivity and ac susceptibility characterization of our sample, (a) Temperature 
dependence of the resistivity of our sample, (b) Temperature dependence of the ac susceptibility 
of our sample (red line - real part X\ blue line - imaginary part X ”). The frequency and amplitude 
of the driving field are respectively 7 Hz and 0.64 Oe. 
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FIG. S2. Summary of time-resolved data. Color-coded maps of differential reflectivity spectra 
as a function of pump-probe delay time at different temperatures and comparable fluences of ~1 
m J / cm 2 . 
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FIG. S3. Summary of spontaneous Raman data. Spontaneous Raman spectra as a function 
of temperature at 2.4 eV (a) and 3.1 eV (b) excitation energies. In the color code, black denotes 
Ty. Red arrows indicate the modes that are nominally optically-active only below Ty. The modes 
are labeled according to the symmetry of their cubic counterparts, (c) Dependence on temperature 
of energy and width of the A g mode in the monoclinic structure that originates from the folding 
of the A 2 mode of the cubic structure from the A- to the T-point. The average width of the mode 
corresponds to a phonon lifetime of 620T50 fs. 
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FIG. S4. Singular value decomposition (SVD) of the differential reflectivity matrix 
and fit of the canonical traces with model functions. Canonical time traces ui(t) (a) and 
U 2 (t) (b) obtained from our SVD algorithm. The fitting functions are comprised of the sum of 
two exponential decays, a step function and two damped coherent oscillations, convolved with a 
Gaussian response function. Time (c) and energy dependence (d) of the physical traces obtained 
from SVD, respectively, {Uj(t)} and (Vj(E)}, which include relaxations (i=1,2,3) and damped 
coherent oscillations (i=4,5). 
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FIG. S5. Comparison between data and calculations of the optical response. The imag¬ 
inary part of the permittivity, € 2 , at 15 K from the ellipsometry measurements of Ref. 4 (open 
circles) is plotted together with the theoretical predictions from our GGA+U calculations for dif¬ 
ferent values of the U parameter and exchange coupling ,1 = 1 eV (lines). 
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FIG. S6. Assignment of the interband transitions below 3 eV. Comparison between the 
absorptive part of the optical conductivity from experiments (open circles, Ref. 4) and our GGA+U 
calculations for U = 4 eV and J = 1 eV (black line). The separate contributions from the interband 
transitions below 3 eV are also plotted (colored lines) and labeled according to the predominant 
character of the initial and final states. 
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FIG. S7. Energy dispersion of the initial and final states of the interband transitions in 
Fig. S6. Density of states and band structure of the minority (a,c,d) and majority spin channel 
(b) of the low-temperature phase of magnetite in the approximate symmetry P2/c obtained from 
our GGA+U calculations with U = 4 eV and J = 1 eV. The initial and final states of the interband 
transitions in Fig. S6 are here highlighted in different colors, namely, peak 1 - Fe^ tig —> Fe^ 
t' 2 g (a), peak 2 - Fe^ + e g -» Fe^ + e (b), peak 3 - Fe^ + t 2 g —> Fe^ + e g (c) and peak 4-02 p , Fe^ + 
3 d —f Fe"^ - t2g (d). 
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